#include "Types.hpp"

__kernel void setCellsToMinus(__global int* cells,
	      		      const uint NCells)
{ 
    uint i = get_global_id(0);
    if (i >= NCells) return;
    
    cells[i]=-1;	
}

__kernel void sortParticleIDs(__global int* particle_ids,
	      		      const uint N)
{
    uint i = get_global_id(0);
    if (i >= N) return;
    
    particle_ids[i]=i;
}

__kernel void assembleCells(const __global real4* positions,
	      		      __global int* part_ids,
			      __global int* cells,
			      const uint N,
			      const uint nx,
			      const uint ny,
			      const uint nz,
			      const real len_x,
			      const real len_y,
			      const real len_z)
{
    uint i = get_global_id(0);
    if (i >= N) return;
    
    int cx = positions[i].x/len_x;
    int cy = positions[i].y/len_y;
    int cz = positions[i].z/len_z;

    int cellID=cx*ny*nz+cy*nz+cz;
    part_ids[i]=atom_xchg(cells+cellID,i) ;
}






__kernel void densityPressureCalculation(const __global real4* positions,
	      			 __global real* densities,
				 __global real* inverseDensities,
				 __global real* pressures,
				 __local real* densitySums,
				 const __global int* cells,
			       	 const __global int* part_ids,
				 const __global real* masses,
	      		       const uint N,
			       const uint nx,
			       const uint ny,
			       const uint nz,
			       const real len_x,
			       const real len_y,
			       const real len_z,
			       const int reflectx,
			       const int reflecty,
			       const int reflectz,
			       const real xmin,
			       const real ymin,
			       const real zmin,
			       const real xmax,
			       const real ymax,
			       const real zmax,
			       const real Cdens,
			       const real k_gas,
			       const real rho0,
			       const real h_smooth2)
{
  uint i = get_global_id(0);
  if (i >= N) return;
  
  int cx = positions[i].x/len_x;
  int cy = positions[i].y/len_y;
  int cz = positions[i].z/len_z;
    real lx = xmax -xmin;
    real ly = ymax -ymin;
    real lz = zmax -zmin;


  //real Cdens = 315.0/(64.0*pow(h_smooth,9)*PI);

  densitySums[get_local_id(0)] =(real)0.0;
  for(int l=-1;l<=1;++l)
  {
	for(int m=-1;m<=1;++m)
	{
		for(int n=-1;n<=1;++n)
		{
			int neighbor_cx=cx+l;
			int neighbor_cy=cy+m;
			int neighbor_cz=cz+n;

			int flagX=0;
			int flagY=0;
			int flagZ=0;
			
			if(reflectx==0)
			{
				if(neighbor_cx<0)
				{
				     neighbor_cx=nx-1;
				     flagX=1;
				}
				else if(neighbor_cx>nx-1)
				{
			     	     neighbor_cx=0;
				     flagX=1;
				}
			}
			else
			{
				if(neighbor_cx<0 || neighbor_cx>nx-1)
					continue;		
			}	
			
			if(reflecty==0)
			{
				if(neighbor_cy<0)
				{
				     neighbor_cy=ny-1;
				     flagY=1;
				}
			     	else if(neighbor_cy>ny-1)
				{
			     	     neighbor_cy=0;
				     flagY=1;
				}
			}
			else
			{
				if(neighbor_cy<0 || neighbor_cy>ny-1)
					continue;   			     	
			}			
			
			if(reflectz==0)
			{
				if(neighbor_cz<0)
				{
				     neighbor_cz=nz-1;
				     flagZ=1;
				}
				else if(neighbor_cz>nz-1)
				{
			    	     neighbor_cz=0;
				     flagZ=1;
				}			
			}
			else
			{
				if(neighbor_cz<0 || neighbor_cz>nz-1)
			        	continue;
			}

			int cellID=neighbor_cx*ny*nz+neighbor_cy*nz+neighbor_cz;

			int j=cells[cellID];
			while(j!=-1)
			{
                                real4 dij = positions[j] - positions[i];
                                //For the periodicity
                                if(flagX==1)
                                   dij.x=dij.x-round(dij.x/lx);
                                if(flagY==1)
                                   dij.y=dij.y-round(dij.y/ly);
                                if(flagZ==1)
                                   dij.z=dij.z-round(dij.z/lz);
                                dij.w=0.0;
				
				real norm = dij.x*dij.x+dij.y*dij.y+dij.z*dij.z;
                                if(norm<=h_smooth2)
				{

                                        real addDen=masses[j]*(h_smooth2-norm)*(h_smooth2-norm)*(h_smooth2-norm);
                                        //if (i!=j)
                                        densitySums[get_local_id(0)] += addDen;
                                }
				j=part_ids[j];
			}
		}
	}
  }
  densities[i] = Cdens*densities[i];// From the Boundary Conditions

  densities[i] += Cdens*densitySums[get_local_id(0)];
  inverseDensities[i] = 1.0/densities[i];
  pressures[i] = k_gas*(densities[i]-rho0);

//printf("%g",pressures[i]);
}

__kernel void accelerationCalculation(const __global real4* positions,
                               __global real4* accelerations,
                               __local real4* forcePressSums,
			       __local real4* forceViscoSums,
			       const __global int* cells,
			       const __global int* part_ids,
			       const __global real4* velocities,
			       const __global real* pressures,
			       const __global real* inverseDensities,
			       const __global real* masses,
			       const real4 gravity,
			       const uint N,
			       const uint nx,
			       const uint ny,
			       const uint nz,
			       const real len_x,
			       const real len_y,
			       const real len_z,
			       const int reflectx,
			       const int reflecty,
			       const int reflectz,
			       const real xmin,
			       const real ymin,
			       const real zmin,
			       const real xmax,
			       const real ymax,
			       const real zmax, 
			       const real Cpress,
			       const real Cvisco,
                               const real h_smooth,
                               __global real* bdflag)
{
  uint i = get_global_id(0);
  if (i >= N) return;
  
  int cx = positions[i].x/len_x;
  int cy = positions[i].y/len_y;
  int cz = positions[i].z/len_z;
  real lx = xmax -xmin;
  real ly = ymax -ymin;
  real lz = zmax -zmin;


  //real Cpress=22.5/(PI*h_smooth*h_smooth*h_smooth*h_smocut_roth*h_smooth*h_smooth);
  //real Cvisco=45.0*dynVisco/(PI*h_smooth*h_smooth*h_smooth*h_smooth*h_smooth*h_smooth);

  forcePressSums[get_local_id(0)] =(real4)0.0;
  forceViscoSums[get_local_id(0)] =(real4)0.0;

 for(int l=-1;l<=1;++l)
  {
	for(int m=-1;m<=1;++m)
	{
		for(int n=-1;n<=1;++n)
		{
			int neighbor_cx=cx+l;
			int neighbor_cy=cy+m;
			int neighbor_cz=cz+n;

			int flagX=0;
			int flagY=0;
			int flagZ=0;
			
			if(reflectx==0)
			{
				if(neighbor_cx<0)
				{
				     neighbor_cx=nx-1;
				     flagX=1;
				}
				else if(neighbor_cx>nx-1)
				{
			     	     neighbor_cx=0;
				     flagX=1;
				}
			}
			else
			{
				if(neighbor_cx<0 || neighbor_cx>nx-1)
					continue;		
			}	
			
			if(reflecty==0)
			{
				if(neighbor_cy<0)
				{
				     neighbor_cy=ny-1;
				     flagY=1;
				}
			     	else if(neighbor_cy>ny-1)
				{
			     	     neighbor_cy=0;
				     flagY=1;
				}
			}
			else
			{
				if(neighbor_cy<0 || neighbor_cy>ny-1)
					continue;   			     	
			}			
			
			if(reflectz==0)
			{
				if(neighbor_cz<0)
				{
				     neighbor_cz=nz-1;
				     flagZ=1;
				}
				else if(neighbor_cz>nz-1)
				{
			    	     neighbor_cz=0;
				     flagZ=1;
				}			
			}
			else
			{
				if(neighbor_cz<0 || neighbor_cz>nz-1)
			        	continue;
			}
			int cellID=neighbor_cx*ny*nz+neighbor_cy*nz+neighbor_cz;
			int j=cells[cellID];
			while(j!=-1)
			{
                                real4 dij = positions[j] - positions[i];
				//For the periodicity
				if(flagX==1)
                                   dij.x=dij.x-round(dij.x/lx);
				if(flagY==1)
                                   dij.y=dij.y-round(dij.y/ly);
				if(flagZ==1)
                                   dij.z=dij.z-round(dij.z/lz);
      				dij.w=0.0;
				
				real norm = sqrt(dij.x*dij.x+dij.y*dij.y+dij.z*dij.z);
                                if(norm<=h_smooth)
				{
					//Pressure calculation
					real4 Fpress = masses[j]*(pressures[i]+pressures[j])*inverseDensities[j]*(h_smooth-norm)*(h_smooth-norm)*(h_smooth-norm)*dij/norm;
					//Viscosity calculation
					real4 Fvisco = masses[j]*(velocities[j]-velocities[i])*inverseDensities[j]*(h_smooth-norm);
					if (i != j)
					{
       				   	   forcePressSums[get_local_id(0)] += Fpress;
					   forceViscoSums[get_local_id(0)] += Fvisco;
					} 								   	
				}
				j=part_ids[j];
			}
		}
	}
  }

if (bdflag[i]==1)
  accelerations[get_global_id(0)] = accelerations[get_global_id(0)]*masses[get_global_id(0)];
else
  accelerations[get_global_id(0)] = (Cvisco*forceViscoSums[get_local_id(0)]-Cpress*forcePressSums[get_local_id(0)])+gravity;


}

real interpolateBoundaryDensity(real riwall ,  __global real* b_density, real cut_r, int Num_B_particles )
{
// riwall is the Actual distance from the Wall. Compute the index where it may lie
    int i= (float)(riwall/cut_r*(Num_B_particles-1));
   // printf("Boundary Position %d %g\n", i, b_density[i]);
    return b_density[i]; // return the avergae of the neighbouring two particles
}
__kernel void setBoundaries( __global real4* positions,
	      		      const uint N,
			      const real xmin,
			      const real xmax,
			      const real ymin,
			      const real ymax,
			      const real zmin,
			      const real zmax,
				  const int reflectx,
			       const int reflecty,
			       const int reflectz,
			       __global real* densities,
                               __global real4* accelerations,
                               const __global real* b_density,
                               const int Num_B_particles,
                               const real cut_r,
                               const real dr,
                               const real dt,
                               __global real* bdflag )
{
    uint i = get_global_id(0);
    if (i >= N) return;
 densities[i]=0;
 accelerations[i].x=0;
 accelerations[i].y=0;
 accelerations[i].z=0;
 bdflag[i] = 0;

        real left_w= positions[i].x-xmin;
        real right_w= xmax - positions[i].x;
        real bottom_w=  positions[i].y-ymin;
        real top_w= ymax - positions[i].y;
        real back_w=  positions[i].z-zmin;
        real front_w= zmax - positions[i].z;


	if (reflectx==0)
	{
                if(left_w<0)
		positions[i].x+=(xmax-xmin);
                else if(right_w<0)
		positions[i].x-=(xmax-xmin);
	}

        else if (reflectx==1)
        {
            if(left_w <=cut_r && left_w >=0 )
            {
                densities[i]+= interpolateBoundaryDensity(left_w,b_density,cut_r,Num_B_particles);
                if (left_w < dr)
                {
                    bdflag[i]=1;
                    accelerations[i].x = (dr - left_w)/(dt*dt);
                }
             }
            else if (right_w<=cut_r && right_w >=0 )
            {
                densities[i]+= interpolateBoundaryDensity(right_w,b_density,cut_r,Num_B_particles);
                if (right_w < dr)
                {
                    bdflag[i]=1;
                    accelerations[i].x = (- dr + right_w)/(dt*dt);
                }
             }

        }

	if (reflecty==0)
	{
    
                if(bottom_w<0)
		positions[i].y+=(ymax-ymin);
                else if(top_w<0)
		positions[i].y-=(ymax-ymin);
	}

        else if (reflecty==1)
        {
            if(bottom_w <=cut_r && bottom_w >0)
            {
                densities[i]+= interpolateBoundaryDensity(bottom_w,b_density,cut_r,Num_B_particles);
                if (bottom_w < dr)
                {
                    bdflag[i]=1;
                    accelerations[i].y = (dr - bottom_w)/(dt*dt);
                }

            }
            else if (top_w<=cut_r && top_w >0)
            {
                densities[i]+= interpolateBoundaryDensity(top_w,b_density,cut_r,Num_B_particles);
                if (top_w < dr)
                {
                    bdflag[i]=1;
                    accelerations[i].y = (-dr +top_w)/(dt*dt);
                }
            }
          }


	if (reflectz==0)
	{
                if(back_w <0)
		positions[i].z+=(zmax-zmin);
                else if(top_w<0)
		positions[i].z-=(zmax-zmin);
	}

        else if (reflectz==1)
        {
            if(back_w <=cut_r && back_w >=0)
            {
                densities[i]+= interpolateBoundaryDensity(back_w,b_density,cut_r,Num_B_particles);
                if (back_w < dr)
                {
                    bdflag[i]=1;
                    accelerations[i].z = (dr - back_w)/(dt*dt);
                }
            }
            else if (front_w<=cut_r && front_w >=0)
            {
                densities[i]+= interpolateBoundaryDensity(front_w,b_density,cut_r,Num_B_particles);
                if (front_w < dr)
                {
                    bdflag[i]=1;
                    accelerations[i].z = (-dr + front_w)/(dt*dt);
                }

            }
        }

}


__kernel void updatePositions( __global real4* positions,
	      		       __global real4* old_forces,
			      const __global real4* forces,
                              const __global real4* velocities,
                              const __global real* inv_2masses,
                              const real dt,
                              const uint N)
{
    uint i = get_global_id(0);
    if (i >= N) return;

    positions[i]=positions[i]+dt*velocities[i]+forces[i]*dt*dt*inv_2masses[i];
    old_forces[i]=forces[i];
}

__kernel void updateVelocities(__global real4* velocities,
	      	               const __global real4* old_forces,
                               const __global real4* forces,
			       const __global real* inv_2masses,
                               const real dt,
                               const uint N)
{
    uint i = get_global_id(0);
    if (i >= N) return;
    velocities[i]=velocities[i]+(old_forces[i]+forces[i])*inv_2masses[i]*dt;
}

__kernel void eulerIntegrator( __global real4* positions,
	      		      __global real4* velocities,
	 		      const __global real4* accelerations,
			      const real dt,
			      const uint N)
{
	const uint pID = get_global_id(0) ;
	if(pID >= N) 
	  return;
		
	positions[pID] = positions[pID] + velocities[pID]*dt;
	velocities[pID] = velocities[pID] + accelerations[pID]*dt;
}

__kernel void yusukeIntegrator( __global real4* positions,
	      		      __global real4* velocities,
	 		      const __global real4* accelerations,
			      const real dt,
			      const uint N)
{
	const uint pID = get_global_id(0) ;
	if(pID >= N) 
	  return;

	positions[pID] = positions[pID] + velocities[pID]*dt +accelerations[pID]*dt*dt*0.5;
	velocities[pID] = velocities[pID] + accelerations[pID]*dt;
}
